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ABSTRACT 

Results of numerical 3D simulations of propagation of acoustic waves inside 
the Sun are presented. A linear 3D code which utilizes realistic OPAL equation 
of state was developed by authors. Modified convectively stable standard so- 
lar model with smoothly joined chromosphere was used as a background model. 
High order dispersion relation preserving numerical scheme was used to calcu- 
late spatial derivatives. The top non-reflecting boundary condition established in 
the chromosphere absorbs waves with frequencies greater than the acoustic cut- 
off frequency which pass to the chromosphere, simulating a realistic situation. 
The acoustic power spectra obtained from the wave field generated by sources 
randomly distributed below the photosphere are in good agreement with obser- 
vations. The influence of the height of the top boundary on results of simulation 
was studied. It was shown that the energy leakage through the acoustic potential 
barrier damps all modes uniformly and does not change the shape of the acoustic 
spectrum. So the height of the top boundary can be used for controlling a damp- 
ing rate without distortion of the acoustic spectrum. The developed simulations 
provide an important tool for testing local helioseismology. 

Subject headings: Sun: oscillations — sunspots 



1. Introduction 

Solar 5-min. oscillations are excited by turbulent convection (downdrafts) in subsurface 
layers of the Sun. These oscillations consist of acoustic and surface gravity waves with 
frequencies in the range of 2^8 mHz and in a wide range of wave numbers. The observed 
oscillations can be used for reconstruction of internal structure of the Sun by methods of 
helioseismology. There are several methods of investigation of the interaction of traveling 
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acoustic waves wit h small perturbations to the backgro und state. One of them is the time- 
distance approach ( Duvall et al. 1993 : Kosoviche vl 1 1 9 96l ) . The key concept of this method is 
measuring and inversion of wave travel times. Propagation of acoustic waves in this approach 
is calculated using variou s approximations to the wave equations, such a.s the ray theory or 



first Born approximat i on (IKosovichev fc Duvall 1119971 : iKosovichev et al. Il2000l : iJensen et al. 



200ll : ICouvidat et al.l 12004 120061). The s e app r oximations hav e been tested using simple 
models for point sources (e.g. iBirch et al.l ( 200ll ): lBirch fc Felderl ( 20041 )). but not for realistic 
solar conditions, e.g. realistic stratification and random excitation sources. Such tests, which 
require direct numerical simulations, are extremely important for validating the inferences 
from time-distance helioseismology and other local helioseismology methods. 

There are two main directions in numerical simulations of solar oscillations and waves. 
The first one is to use realistic non-linear simulations of solar convection. In such simula- 
tions, the waves are naturally excited b y convective motio ns. These simulations reproduce 
quite well the solar oscillation spectrum IStein et al.l (12004) . and have been used for testing 
time-distance helioseismology Georgobiani et al'] ~( 2006 ). The second approach is based on 
linearized Euler equations describing wave propagation for a given background state. The 
background state can be taken from non-linear numerical simulations or by perturbing the 
standard solar model. In this paper, we describe a numerical method and initial simulation 
results developed for the second approach. 

We developed a 3D code which utilizes a realistic physics and accurately simulates re- 
flection of acoustic waves from the top boundary. A realistic equation of state which takes 
into account partial ionizations and various corrections is used. The equation of state was 
calculated by interpolation of the OPAL tables. In Section 2.1 we give a detailed description 
of the underlying physics. The main attention is paid to developing a consistent procedure 
for obtaining a convectively stable background model and establishing realistic top boundary 
conditions based on the Perfectly Match Layer (PML) method. In Section 2.2 we describe 
a semi-discrete numerical scheme of high order which preserves dispersion relations of the 
continuous problem. The main attention is paid to developing stable high-order numerical 
boundary conditions consistent with the finite-difference scheme, in which the dispersion re- 
lation is preserved for the inner mesh points of the computational domain. In Section 3, we 
compare the numerical and analytical solutions of various ID test problems with an isother- 
mal background model to validate the code, investigate accuracy of the numerical scheme, 
and the non-reflecting boundary conditions in a gravitationally stratified medium. In Section 
4, we present results of numerical three dimensional simulations of acoustic wave field for 
a standard solar background model. We used different types of single or multiple acoustic 
sources: z-component of force or pressure, point or distributed with different time depen- 
dence. The main goals are to study properties of solar waves for various excitation sources 
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and generate artificial wave fields for testing accuracy of the Born and ray approximations 
and local helioseismic diagnostics of the solar interior, currently used for SOHO/MDI and 
GONG data. The results of this testing will be presenting in future papers. The numerical 
simulations are carried out on parallel supercomputers at NASA Ames Research center. 



2. Code description 
2.1. Physical background 



The propagation of adiabatic acoustic waves below the solar photosphere is described 
by the following system of linearized Euler equations: 
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where u', v', w' are the perturbations of x, y, z velocity components, p' and p' are the den- 
sity and pressure perturbations correspondingly. Quantities with subscript such as the 
pressure po, the density po, and the gravitational acceleration go correspond to the back- 
ground r eference model and depend only on radius. Christensen-Dalsgaard's standard solar 
model S flChristensen-Dalsgaard et al.lll996l ) with smoothly joined chromosphere provided by 



Vernazza et al 



fll976l ) was chosen as such a model. To close the system we use an adiabatic 



relation between Eulerian variations of pressure p' and density p' 
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where Oq = Fipo/po is the square of sound speed, Fi = {dlogp/ d log p) ad is the adiabatic 
exponent, Nq is the Brunt-Vaisala frequency, C,z is the vertical displacement. So, from the 
background model we need only profiles of the sound speed and Brunt-Vaisala frequency. 

The background state is convectively unstable, especially just below the photosphere, 
where the temperature gradient is super-adiabatic and convective motions are very intense 
and turbulent. This leads to the instability of the solution of linear system ([1]). Convection 
instability is developed on a time scale of 30-40 minutes of solar time. Simulation of the 
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acoustic spectrum of solar oscillations needs to perform calculations on a time interval of 
5-8 hours of solar time, and this instability will badly distort the result. To make the 
background model stable against convection we slightly modified profiles of pressure and 
density in a thin (500 km in depth) super-adiabatic layer just below the photosphere. The 
condition for stability against convection requires that the square of Brunt- Vaisala frequency 
has to be positive 

N\r)=g(^^-^]>0, (3) 



Fi dr dr 

where r is the distance from the center of the Sun. The profile of A^^ near the solar surface 
is shown in the bottom left pane of Figure [1] by the solid line. If we replace negative values 
by zero (or small positive) value we guarantee a stability of the modified model against 
convection. Now we can recalculate the pressure and density from the modified profile of 
^mod- Combining equation ([3]) with the condition of hydrostatic equilibrium, we get the 
following boundary value problem for p and p: 

1^ = 3 cc^Ld 

pdz (? (yf ' 

dp (4) 
dl = -^^' 

0<z<L, p(0) = po(0), p{L)=p^{L), p{L)=po{L), 

where L is the depth of the computational domain, z is the vertical coordinate which is 
counted off from the bottom of the domain. We introduced a free parameter a that has 
to be determined, to make the modified model closer to the original one. We established 
boundary conditions for the system (jlj) setting the density at the top and bottom boundaries 
equal to solar values. Parameter a does not change the condition of convective stability if 
it remains positive. Introducing of an additional parameter permits us to establish yet 
another boundary condition and fix the pressure at the top boundary. So the recipe to 
build a convectively stable background model close to the standard one is the following. We 
smoothly join the density profiles of the standard solar model and the chromosphere, obtain 
the pressure profile from the condition of hydrostatic equilibrium, calculate square of the 
modified Brunt- Vaisala frequency profile -/Vmod; replacing negative values by a zero or small 
positive value, and substitute it into the right hand side of (jlj). Parameter a, profiles of 
density and pressure of the modified convectively stable model are obtained as a solution of 
the eigenvalue problem (jlj). The profile of Fi , needed for the sound speed profile, is found 
from the realistic OPAL equation of state fjRogers et al.jjl996j ). Vertical profiles of pressure 
Po, density po, sound speed oq, square of Brunt- Vaisala frequency 7V^, adiabatic exponent 
Fi, and acoustic cut-off frequency 

.2 «o A o^^^ iT-i dlogpo 
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for both models (standard with joined chromosphere and modified convectively stable) are 
shown in Figure [H The depth of the domain from the photosphere equals 30 Mm, the 
height of the chromosphere equals 2 Mm, a = 0.861. The dashed curves represent profiles 
in the convectively stable model, the solid lines show profiles of the standard solar model. 
The thin vertical line marks the position of the fitting point of the chromosphere and the 
standard solar model. The background model, obtained with the procedure described above, 
is convectively stable and self-consistent. 



where q = {p' , pou' , pov' , pow')^ is the vector of independent variables, S{q,t) is the source 
term containing the gravity term and acoustic sources which depend on time explicitly. The 
source term does not contain spatial derivatives. Explicit expressions for components of 
vectors F, G, and H can be easily found from the system ([T]). We used a semi-discrete 
numerical scheme. In semi-discrete approach the space and time discretization processes 
are separated. First the spatial discretization is performed, leaving the problem continuous 
in time. The spatial derivatives have been approximated by the finite difference method 
reducing the system of partial differential equations to the system of ordinary differential 
equations 



2.2. Numerical algorithm 



The system ([T]) is written in a conservative form 
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which can be solved by any stable time advancin g method . We used four stage, 3rd order 
strong-stability-preserving Runge-Kutta method (jShull2002l ) with Courant number c = 2: 
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f^ + ^AtL{q('\r + Ut). 



The source term S{q, t) which depends on time explicitly does not require a special treat- 
ment. 



High-order dispersion-relation-preserving (DRP) scheme developed by iTam fc Webb 



( 119931 ) was used for spatial discretization. Coefficients ai of the finite difference scheme 
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are chosen from the requirement that the difference between Fourier transform of the nu- 
merical scheme and Fourier transform of the spatial derivative has to be minimal. Taking 
the Fourier transform of both sides of iQ one can get the effective wave number k^ff of the 
Fourier transform of the numerical scheme (19]) 
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An assumption that the integral error E of Fourier transform of the finite difference scheme 
([9]) is minimal for waves with wavelength A > 4Ax leads to the following equation 
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are obtained from a requirement that the numerical scheme ([9]) approximates a spatial deriva- 
tive with the 4th order. Now calculation of the coefficients of 7-dot symmetrical stencil is 
straightforward. The explicit expressions for the coefficients are the following: 
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The plots of numerical wave number k^ffAx versus kAx for different numerical schemes are 
shown in Figure [21 Dotted, dash-dotted, dashed, and solid curves represent classic 2nd, 4th, 
6th, and DRP 4th order schemes correspondingly. One can see that the 4th-order DRP 
scheme describes short waves more accurately than the classic 6th-order scheme. 

Waves with wavelength less than 4 Ax are not resolved by the numerical scheme 
They lead to point-to-point oscillations of the solution that can cause a numerical instability. 
Such waves have to be filtered out. We used the following digital filter of the 6th-order to 
eliminate unresolved short wave component from the solution: 



fsmix) = f{x) - crfD{x) = f{x) - CT/ X dmf{x + mAx), 



(15) 



where / is the original grid function, fsm is the filtered grid function, D{x) is the damping 
function, o"/ is the constant between and 1, determining the filter strength. The frequency 
response function G{k) of the filter relates the Fourier images of the original / and filtered 
fsm grid functions as follows fsmik) = G{k)f{k). In this paper, the coefficients of the 
digital filter have been chosen in such a way that 
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Coefficients dm of the digital filter are symmetric: 
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Using the technique proposed by ICarpenter et al.l fll993l ) we have found a stable 3rd-order 
boundary closure of the explicit dispersion-relation-preserving inner scheme (Q with coeffi- 
cients given by f|T^ . Obtained boundary closure has summation- by-parts properties. This 
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approach is based on the imphcit Pade approximation of spatial derivatives near the bound- 
aries 

= (18) 
where matrices P and Q satisfy the following conditions: 

1. P is symmetric non-singular matrix (P = P"^), 

2. P is positive-definite matrix {U^ PU > for V C/), 

3. Q is almost skew-symmetric matrix, except corner top left and bottom right elements 
(Q + Q^ = |go,o|diag(-l,0,...,l)) 

4. qN,N > 0, go,o = -qN,N- 



Taking into account these properties, one can write explicitly the top left corners of matrices 
P and Q: 
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where coefficients a^ of the inner scheme are defined in (fT4l) . Expanding the left and right- 
hand sides of (1181) in Taylor series at the top boundary and equating terms of the same order 
of Ax, one can obtain a system of linear equations for coefficients Pij and qij. Not all of these 
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equations are independent, so the solution depends on two free parameters pss and P23: 

83 34016 - 149737r 
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To satisfy a condition of positive definiteness it is sufficient to clioose matrix elements P33 and 
P23 in sucli a way, tliat tfie signs of coefficients of a characteristic polynomial ahernate. How- 
ever, this property does not guarantee a boundedness of the solution for all times, which is 
called asymptotic stability. To make a solution be bounded for all times, all eigenvalues of the 
spatial discretization operator Ln^j from Eq.(j7]), incorporated with the bou ndary conditions. 



must have non-positive real parts. Details of this procedure can be found in (jCarpenter et al. 



I993I ). Due to complexity of the original 3D problem, we have tested a stability of the scheme 
on ID advection problem. Distribution of eigenvalues of the DRP spatial discretization op- 
erator in the complex plane for different choices of the pairs of coefficients (p23! P33) is shown 
in Fig. O Symbols plus correspond to the scheme P23 = 1/30, P33 = 31/32, which does 
not exhibit asymptotic stability. Circles and crosses represent choices (1/80, 125/128) and 
(-1/10, 65/64) of coefficients (^23,^33) correspondingly. Both schemes are asymptotically 
stable. 

Acoustic sources have been added to the right-hand side of equations ([1]). We used 
sources of two types. If we add a scalar function ^{x,y, z,t) to the right-hand side of z- 
component of momentum equation, this term can be combined with the gravity term and 
interpreted as a source of z-component of force. If we add a gradient of a scalar func- 
tion V$(x, z,t) to the right-hand side of all momentum equations, these terms can be 
combined with components of the pressure gradient and interpreted as a pressure source. 
Acoustic sources are spatially localized and have finite lifetime. Spatial dependence is given 
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by Gaussian spherically symmetric function with semi-width of 2-3 grid nodes. We ex- 
perimented with two different time dependencies of acoustic sources: one period of sin 
function sin[u;(t — to)], to < t < + I'k joj and Ricker's wavelet (1 — 2x^)e~^^, x = 
[uj{t — to)/2 — tt], tQ < t < to + in/uj. Such time dependencies were chosen because 
such sources are not monochromatic and have spectral power localized around central fre- 
quency uj/271, but spectral power is not too spread out. Single or multiple acoustic sources 
can be added to the right hand side of equations. Multiple sources are randomly distributed 
at some depth (in our case it was 350 km) and initiated independently at arbitrary moments 
of time. Amplitudes and frequencies are randomly distributed on intervals [0, 1] and [2 mHz, 
8 niHz] correspondingly. Thickness of both (top and bottom) PMLs equals 5 grid nodes. 

Calculation of the acoustic spectrum requires long term simulations, so besides the 
asymptotic stability we have to prevent a spurious reflection of acoustic wa v es fro m the 



boundaries back to the computational domain. In this paper we follow iHul (119961 ). who 
proposed a procedure to construct the PML for the Euler equations. It can be proofed that 
for a homogeneous medium and uniform mean flow without gravity PML absorbs waves 
without reflection for any angle of incidence and frequency. We set non-refiecting boundary 
conditions based on the PML at the top and bottom boundaries of the domain. The lateral 
boundary conditions are periodic. Inside the PML independent variables q are split into 
sub- components qi,q2, q^ such that q = qi + q2 + Qs- Thus, in the PML 3D system ([1]) is 
split into ID+ID+ID system of coupled locally one dimensional equations 

9Qi I dFjq) _ 

dt dx 

%+^ = 0. (21) 

at oy 

dqs dH{q) 

-g^ + ^^ = S{q,t)-a.q„ 

where Ata^ = 0.05 + amax{Z / D)"^ is the damping factor, Z is the vertical coordinate inside 
the PML counted off from the interface of the PML with the inner region, D is the depth of the 
PML. Values of amax at the top and bottom boundaries are 0.3 and 1.0 correspondingly. In 
the paper of F.Q. Hu a quadratic dependence of a on the coordinate Z is used. We were forced 
to add a small constant term to stabilize the PML in the presence of gravity. It is important 
to note, that vectors F, G, and H depend only on unsplit variable q. Although q^, q2, and 
q^ are not defined outside the PML, the variable q, which is used for calculation of the spatial 
derivatives, is defined everywhere in the computational domain. Hence, inside the PML near 
the interface with the inner region we can use the same centered stencil as for the inner points. 
Near the top and bottom boundaries the implicit Pade approximation ( fT8|) is used which 
guarantees numerical stability of the scheme. We smoothly joined the chromospheric model 
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provided by IVernazza et al.l (119761 ) with the top of the standard solar model by Christensen 
Dalsgaard and established the top non-reflecting boundary condition based on the PML 
in the chromosphere above the temperature minimum. This simulates a realistic situation 
when not all waves are reflected by the photosphere. Waves with frequencies higher than the 
acoustic cut-off frequency pass through the photosphere and will be absorbed by the PML 
layer. 



3. Numerical examples 

For validation of the code we chose ID initial boundary value problem (IBVP) for 
linearized Euler equations with constant gravity Qq = const: 

< X < 1, t > 0, 
p'(0,t) = p'(l,t)=0, 

p'(x,0) = /(x), po{x)u'{x,0)=0, (22) 

X 

po{x)^{x, 0) = - j f{ri) dr], 



where ^ is the displacement. We need to know C, to calculate the Eulerian perturbation of 
pressure p'. Actually, the combination p^^ is used as a variable, so we group them together in 
equations. Initial conditions for ^ must be consistent with the initial conditions for p'. They 
are related by the continuity equation. Waves are adiabatic, the background model po,Po 
is hydrostatic and isothermal Po/Po = const. The last equation in the system (12^ is the 
adiabatic relation ([2]) written for the isothermal background model. Variable x represents 
here the depth from the surface. The system fl22|) is written in the same conservative form 
as the original system ([1]). This problem was chosen for testing the code because it shows all 
characteristic behavior of the realistic solution and yet not too complicated and can be solved 
analytically. Formally, the variable po^ can be eliminated using the continuity equation in 
ID, and system fl22]) can be reduced to the system of two equations. However, in 3D p^C, 
cannot be eliminated. To have the test example as close to the real possible, we 

left this variable and solved numerically the full system (122|) . Analytical solution of these 
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equations can be obtained by the method of separation of variables 

oo 

t) = e^l"^^ An sin Trnx cos A„aot, 

n=l 
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^(x, t) = e~^^'^^ i?„(sin irnx — 2TinH cos Trnx) cos A„aot 
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An = 2 [ /(r/)e-''/2^^sin7rnr/dr/, 5„ = 



n=l (23) 

2//A„ 



1 + Air^n'^H^ ' 

19o 



The following distribution of density perturbation f{x) was chosen as the initial condition 
for p': 

10^[(x - 0.5)2 - 0.001]2 ifo.4<x<0.6 



fix) = I ' ~ ~ ' ^^^^ 

Solution of the initial boundary value problem (l22l) for different moments of time with 
the initial density distribution given by fl24l) and parameters = 1, 'j = 5/3, Qq = 10, 
At = 2 ■ 10""^, = 200 (number of grid nodes) is shown in Figure HI The left column 
represents the density perturbation, the right one shows the vertical displacement. The 
solid curve is the exact solution fl23|) . The dashed curve represents the low-order (classic 
2/1) numerical solution which uses the 2nd-order classic central difference approximation of 
spatial derivatives for inner points with the one sided Ist-order scheme at the boundaries. 
The high-order numerical solution is indistinguishable from the exact one. It uses the DRP 
4/3 scheme (dispersion-relation-preserving spatial discretization of the 4th order for inner 
points with the stable boundary closure of the 3rd order consistent with the inner scheme). 
The bottom panels give the profiles of density and displacement after reflection from the 
bottom boundary. One can see, that the second order solution approximates the exact one 
well enough before the wave hits the boundary. After this the accuracy of the solution 
switches from the second to the first order. Solution becomes too dispersive which causes 
nonphysical oscillations. The high-order solution based on DRP 4/3 scheme reproduces the 
exact solution well even after 30000-^40000 iterations and 20-^30 reflections from boundaries. 
This test shows that the high-order DRP numerical scheme does not introduce a noticeable 
damping or dispersion even on big intervals of integration. These simulations also test an 
accuracy and stability of the numerical boundary conditions. 

To test the efficiency of the PML for non-uniform isothermal background model we 
compared the numerical solution of fl22p with the PML established at the top boundary with 
the exact solution of the same problem for infinite interval — cxo < x < cxo: 

p'{x,t) = l/(x + aot)e-"°*/2^ + i/(a;-aot)e'^«*/2//_ 
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M^a^t^-{x^/2H) 



(25) 



Bottom boundary condition for the numerical solution remains reflecting, because the bottom 
PML is inconsistent with the initial conditions for ^{x,t). At the bottom boundary in the 
initial moment of time ^(1,0) < (see the top right panel of Figure H]). If we established 
PML at the bottom, it would damp ^ to zero value, generating non-physical perturbations 
near the bottom boundary which corrupt the solution. The analytical solution ( 1251) does 
not contain reflected waves, because all initial perturbations propagate to inflnity. This 
solution can be used as a reference solution for determining the damping properties of the 
top PML. Results are shown in Figure [5] at the moments of time t = 0, 0.2, 0.4, and 0.64. 
The value p' / ^/po (density perturbation with removed exponential factor) is plotted. The 
solid line represents the exact solution ( l25i) . the dash-dotted line represents the numerical 
solution with PML at the top boundary, and the the dashed line represents the exact solution 
(!23|) for the reflecting top boundary. The solid vertical line marks position of the interface 
between the top PML and inner region. The dashed vertical line shows position of the initial 
perturbation. The top PML reduces the amplitude of reflected wave by factor 20-i-40. 

Our original 3D system contains acoustic sources explicitly depending on time in the 
right hand side of the momentum equations. We have tested the code in presence of sources 
on the same problem fl22|) with the pressure source term 
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(poO = Pou', 



p'(x,0) = 0, po{x)u'{x,0) = 0, 



P' = alp' + (7 



l)^o(PoO> 



Po(x)e(x, 0) = 0. 



For test purposes we chose gaussian shaped harmonic source function as follows 
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The system (1261) can be solved analytically by quadratures 

t 1 








G{x, r/, r) = 2 e^^-^^/^H 




sinTrna; sm nnr], 



(28) 



The results of numerical simulations with the source function given by ( 1271) and parameters 
= 120, At = 2 - 10-3, ao = 1, 7 = 5/3, go = 10, hsrc = 0.4, uq = lOvr, a = 0.0178 are 
shown in Figure [61 The non-refiecting boundary conditions are established on the top and 
bottom boundaries for numerical solution. The solid curve represents the exact solution of 
([6]) with zero boundary conditions for p' established at x = and x = 1. The dashed line 
represents DRP 4/3 numerical solution of ([6]) with non-refiecting top and bottom boundaries. 
The vertical dashed line marks the position of the source. The vertical solid line shows the 
position of the interface between the inner region of the computational domain and the non- 
refiecting PML. The numerical solution reproduces the exact one well in the inner region, 
and is effectively damped by the absorbing layer, preventing unwanted reflection from the 
bottom boundary. 

Numerical simulations of propagation of waves in 3D from a single source inside the 
Sun are shown in Figure [71 The Brunt-Vaisala frequency of the standard solar model 
with smoothly joined chromosphere was modified near the surface to make the model sta- 
ble against convection. Such modified model was chosen as the background model. Non- 
refiecting boundary conditions (PMLs) were established at the top and bottom boundaries. 
The top layer was established at the height of 500 km above the photosphere in the region 
of the temperature minimum. This layer absorbs all waves with frequencies higher than 
the acoustic cut-off frequency which pass to the chromosphere and do not affect reflection 
of waves with lower frequencies, because these waves are reflected from layers below the 
photosphere. Lateral boundary conditions are periodic. The computational domain of size 
120x120x50 Mm^ was covered by the uniform grid of size 720x720x300 nodes with spatial 
intervals Ax = Ay = Az = 170 km. The time step At = 1 sec was chosen from stability 
condition. The Gaussian spherically symmetric pulse source of z-component of force 



with a = 0.4 Mm was placed at the depth of 3.4 Mm below the photosphere. For such a 
choice of a semi-width of the source equals approximately 4 grid nodes. The time dependent 




(29) 



- 15 - 



part is just one period of sin-function with ujq = 2.5 mHz. Figure [7] shows snapshots of the 
density perturbation from such a source at t = 11.7 min (left column) and t = 21.7 min 
(right column). The top row represents the vertical slices of the computational domain, the 
bottom row shows the horizontal slices at a height of 350 km above the photosphere. The thin 
horizontal line at ^ = represents the photosphere. The left column shows the disturbance 
from the direct wave, generated by the source. The right column shows the wave reflected 
from the photosphere. The reflected wave front is broader and has less amplitude, because 
our source has finite lifetime (one period of sine) and generates high frequency waves which 
pass through the photosphere. Such waves are absorbed by the top non-refiecting layer and 
do not make a contribution to the amplitude of the reflected wave. 



4. Results and Discussion 



The developed numerical method and code have been used for simulation of the acoustic 
wave field generated by multiple acoustic sources inside the Sun. We found that the height 
of the PML affects absorbing properties of the top boundary and the shape of the acoustic 
spectrum {k-u diagram). Reflection from the top boundary is a wave process. The waves 
are reflected not from the fixed level but from some vertical region. Region with the acoustic 
cut-off frequency greater than the wave frequency > uj acts as a potential barrier for 
such waves. Even if the wave frequency is less than the acoustic cut-off frequency waves 
penetrate to this region with exponentially decaying amplitude. If the thickness of the 
barrier is finite waves can leak through it. This process is similar to the tunneling effect 
in quantum mechanics. This happens in the real Sun as well. We studied behavior of the 
solution for different heights of the top boundary. All depths and heights are calculated from 
the level of the photosphere (r = Rsun in the Christensen-Dalsgaard's standard solar model 
S). The background model varies fast in the region above the temperature minimum. To 
be able to simulate propagation of acoustic waves in the chromosphere we were forced to 
reduce the vertical spatial step to Az = 50 km to preserve the numerical stability. To keep the 
horizontal size of the domain as in previous simulations without significant increasing number 
of grid nodes the horizontal spatial steps were chosen 3 times bigger Ax = Ay = 3Az = 150 
km. To satisfy the Courant stability condition for the explicit scheme, the time step was 
reduced to At = 0.68 sec. The computational domain of size 122.2 Mm x 122.2 Mm x 
32 Mm is covered by the uniform grid of size 816x816x640. Sources of z-component of 
force with random frequencies were randomly distributed at the depth of 350 km. Sources 
are initiated at random moments of time (one source per time step) and depend on time 
as Ricker's wavelet with central frequency from range 2^8 mHz. In the Sun we have some 
damping due to the turbulent viscosity. We simulated this additional viscosity by adding 
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a damping term —adQz to the right hand side of z-momentum equation in the region above 
the photosphere and smoothly fading to zero below. The time dependence of RMS (root 
mean square) wave amplitude averaged along the horizontal plane at the height of 300 km 
above the photosphere for different heights htop of the top boundary and different values of 
the damping coefficient ad is shown in Figure [H The RMS amplitude for the high PML 
established at the height of htop = 1750 km without additional damping cr^ = is shown 
by the solid curve I. The RMS amplitudes for the same height of the top PML and 
= 0.3, 0.6, and 1.0 are plotted by solid lines II, III, and V correspondingly. In the last 
case the RMS amplitude reaches an equilibrium state. The curve IV corresponds to the 
low PML, established at the height of 500 km above the photosphere without additional 
damping in the inner region ad = 0. RMS amplitude reaches an equilibrium state in this 
case as well, because the acoustic modes leak through the acoustic potential barrier and 
their exponential tails reach the top absorbing boundary, which adds an additional damping 
and stabilizes the amplitude. The top boundary with the height of 1750 km is set high 
enough and does not affect modes with frequencies less than the acoustic cut-off frequency. 
Energy is continuously pumped to the system by acoustic sources. The bottom boundary 
is deep enough that some modes resolved by numerical scheme have turning points above 
the bottom boundary. Such modes are trapped in the domain, the total energy increases, 
and the RMS amplitude does not reach an equilibrium state. This distorts the acoustic 
power spectrum and changes the amplitude ratio of trapped modes and modes that can be 
absorbed at the top and/or bottom boundaries. The left panes in Figure [9] show the acoustic 
power spectra obtained from observations (top), simulations without damping with low (500 
km) top boundary (middle), and high (1750 km) top boundary (bottom). The right panes 
show the vertical cuts of corresponding k-u diagrams at / = 584. The bottom left pane 
(high top boundary) shows presence of g- modes in simulations. They appear because our 
background model is convectively stable in the thin layer below the photosphere. In the real 
Sun, this layer is convectively unstable, and thus the g- modes do not propagate. Energy 
leakage through the acoustic potential barrier in the case of low top boundary (middle row) 
damps all modes uniformly and does not change the shape of the acoustic spectrum. So the 
height of the top boundary can be used for controlling the damping rate without distortion 
of the acoustic spectrum. 



5. Conclusion 

Developed linear 3D code for propagation of acoustic waves inside the Sun uses the 
realistic equation of state and realistic non-reflecting boundary conditions which permits to 
simulate accurately reflection of waves from the top boundary. Waves with frequencies less 
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then the acoustic cut-off frequency are reflected from the photosphere, and waves with higher 
frequencies pass to the chromosphere. The top non-reflecting boundary absorbs such waves. 
Estabhshing the top boundary high enough in the chromosphere leads to the trapping of 
some modes in the domain and increasing their amphtudes with time, which distorts the 
shape of the acoustic spectrum. Energy leakage through the acoustic potential barrier in 
case of the low (500 km) top boundary leads to an additional uniform damping of all modes 
which stabilizes the amplitude and does not distort the spectrum. The height of the non- 
reflecting top boundary can be used as a parameter for controlling of the damping rate in 
the system. The acoustic spectrum obtained from simulated wave field shows existence of 
p-, and f-modes. The simulated acoustic spectrum is good agreement with observations. 



This code has been used by (jParchevsky fc Kosovichev 1 12006| ) to model the effects of 



non-uniform spatial distribution of acoustic sources in sunspot regions. Their results showed 
that this effect can explain at least a half of the observed amplitude reduction in sunspots. 
The code can be used for studying details of interaction of waves with inhomogeneities of solar 
structure and for producing artificial data for testing an accuracy of helioseismic inversion 
as well, as for studying of propagation of acoustic waves in the chromosphere and reflecting 
properties of the photosphere. Future simulations will include subsurface flows and magnetic 
field. 
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Fig. 1. — Vertical profiles of the density, pressure, sound speed, adiabatic exponent, Brunt- 
Vaisala frequency, and acoustic cut-off frequency. Solid curves represent profiles for the 
standard solar model with smoothly joined chromosphere, dashed ones show the profiles of 
the convectively stable modified model. 
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Fig. 2. — Effective wave number k^ffAx versus kAx for different numerical schemes. Dotted, 
dash-dotted, dashed, and sohd curves represent classic 2nd, 4th, 6th, and DRP 4th order 
schemes correspondingly. 
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Fig. 3. — Eigenvalues of the DRP spatial discretization operator for scalar advection equation 
on complex plane for different choices of the coefficient p23- Symbols plus correspond to the 
scheme P23 = 1/30, P33 = 31/32 which does not exhibit an asymptotic stability. Circles and 
crosses represent choices of P23 = 1/80, —1/10 and P33 = 125/128, 65/64 correspondingly. 
Both schemes are asymptotically stable. 




Fig. 4. — Solution of the IBVP (l22l) for the isothermal hydrostatic background model. 
Density and displacement are shown on the left and right panes correspondingly. Solid curve 
represents the exact solution, dashed one shows the classic 2/1 numerical solution. The high 
order DRP 4/3 numerical solution is indistinguishable from the exact one. After hitting the 
boundary classic 2/1 solution changes the global order of accuracy to 1, because boundary 
conditions are realized only with the 1st order of accuracy. 
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Fig. 5. — Density perturbation with removed exponential factor p'/ is plotted. The solid 
line represents the exact solution for infinite interval (!25|) . the dash-dotted line represents the 
numerical solution with the PML at the top boundary, and the the dashed line represents 
the exact solution (j23l) for reflecting left boundary. The vertical solid and dashed lines mark 
positions of the PML interface and initial perturbation correspondingly. 
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Fig. 6. — Solution of the IB VP ( l26l) for the isothermal hydrostatic background model in 
presence of the source. The solid curve represents the exact solution with zero boundary 
conditions for p' . The dashed line represents DRP 4/3 numerical solution of (!26|) with 
non-reflecting top and bottom boundaries. Numerical solution is damped effectively by the 
absorbing layer. 
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Fig. 7. — Snapshots of the density perturbation from spherically symmetric Gaussian pulse 
of a z-component of force at t = 11.7 min (left column) and t = 21.7 min (right column). 
The top row represents the vertical slices of the computational domain, the bottom row 
shows the horizontal slices at a height of 350 km above the photosphere. Thin horizontal 
line at z = represents the photosphere. The depth of the source equals 3.5 Mm below the 
photosphere. 
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Fig. 8. — Wave amplitude averaged along the horizontal plane at the height of 300 km 
above the photosphere for different heights of the top boundary and different damping co- 
efficients. The curve I corresponds to the high top boundary, established at 1750 km above 
the photosphere without any additional damping. The curves II, III, V correspond to the 
same boundary conditions, but different values of damping coefficient aa = 0.3, 0.6, 1.0. 
The curve IV corresponds to the low top boundary, established at 500 km without artificial 
damping. 
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Fig. 9. — Acoustic spectra obtained from observations and simulations with different heights 
of the top boundary (the left panes). Starting from the top: observations, simulations with 
htop=500 km, simulations with htop=^750 km. The thin white curves on the left panes show 
position of observational ridges for /, pi, and p2 modes. The right panes show cuts of 
k-uj diagrams from left panes at 1=584. For simulations with high top boundary without 
additional damping (bottom row) acoustic modes trapped in the domain distort the shape 
of the acoustic spectrum. 



